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ABSTRACT 

Published maps of red giant stars in the halo region of M31 exhibit a giant steUar 
stream to the south of this galaxy, as well as a giant "shelf" to the northeast of M31's 
center. Using these maps, we find that there is a fainter shelf of comparable size on 
the western side as well. By choosing appropriate structural and orbital parameters 
for an accreting dwarf satellite within the accurate M31 potential model of Geehan 
et al. (2006), we produce a very similar structure in an A^-body simulation. In this 
scenario, the tidal stream produced at pericenter of the satellite's orbit matches the 
observed southern stream, while the forward continuation of this tidal stream makes up 
two orbital loops, broadened into fan-like structures by successive pericentric passages; 
these loops correspond to the NE and W shelves. The tidal debris from the satellite 
also reproduces a previously-observed "stream" of counterrotating PNe and a related 
stream seen in red giant stars. The debris pattern in our simulation resembles the shell 
systems detected around many elliptical galaxies, though this is the first identification 
of a shell system in a spiral galaxy and the first in any galaxy close enough to allow 
measurements of stellar velocities and relative distances. We discuss the physics of 
these partial shells, highlighting the role played by spatial and velocity caustics in 
the observations. We show that kinematic surveys of the tidal debris will provide a 
sensitive measurement of M31's halo potential, while quantifying the surface density 
of debris in the shelves will let us reconstruct the original mass and time of disruption 
of the progenitor satellite. 

Key words: galaxies: M31 - galaxies: interactions - galaxies: kinematics and dy- 
namics 



1 INTRODUCTION 

Our neighboring disk galaxy M31 has turned out to be 
a well-equipped laboratory for research into the accre- 
tion and disruption of substructure in galactic halos. A 
stellar stream of length ~ 120 kpc is observed to flow 
from its tip south and behind M31 into M31's center 
(Ibata ct al. 2001: McConnach ie et al.l 120031: llbata et all 
[2004; Guhathakurta ot al. 2006.), while many other irregular 
morphological and kinematic features i ndicative of accretion 
are seen in the galac t ic disk and hal ojlFerguson et al ] |2002l: 
iMerrett etHI 120031: llbata e t al. 201^. The giant south- 
ern stream is thought to be formed at the last pericen- 
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hereafter Paper II). This stream holds out the promise of 
the first precise measurements of the gravitational potential 
in M31's halo. This is important for many reasons, includ- 
ing breakingthe disk-halo degeneracy in the mass distribu- 
tion (PaDerJ) and thereby understanding the response of 
the dark halo to the central baryons. However, there ar e 
two obstacles to these measurements at present JPaoer ij) . 
First, the current observational errors on the distance allow 
too much flexibility in the orbits to constrain the potential 
effectively. Second, the stream does not follow a single orbit, 
but has a gradient of energy along the length of the stream, 
"tilting" the stream in phase space relative to the orbit and 
biasing the measurement of the potential. The amount of 
this tilt depends on the location of the progenitor, which 
is currently unknown. Locating the progenitor would both 
specify the stream-orbit tilt and constrain the allowed or- 
bital trajectories, removing both obstacles. 

The "Northeast Shelf", a faint diffuse feature seen to 
the NE side of M31 in the red giant branch (RGB) star- 
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count maps of lFerguson et aO i2002h . has been suggested as 
a possibility for the current location of the prog enitor, or at 
least as the forward co ntinuation of the stream jlbata et alJ 
l2004l:lFont et alJl2006 l). Below we find that there is a simi- 
lar, even fainter shelf or plateau on the other side of M31, 
which we call the "Western shelf". In this paper, we fo- 
cus on the possibility that the forward continuation of the 
stream makes up these NE and W shelves. In this scenario, 
the shelves are a similar phenomenon to the "shell" sys- 
tems seen around many elliptical galaxies. These systems 
are believed to be produced by tidal disruption of a satellite 



Table 1. Position of shelf edges (dotted lines in Figure [Tl. 



:alaxy on a nearly radial orbit ( e.g.. ISchweizeilll98ol: ^uinn] 
1983: iHernauist fc QuinnlllQS^ : iBarnes fc HernauislillQm 



exactly the phenomenon believed to be responsible for the 
creation of M31's southern stream, so it is not unexpected 
that shell-like features should be present. These coherent 
structures are rich in information. For example, they allow 
sensitive tests of the potential of the host galaxy without the 
need for distance information (Mcrrificld & Kuiikcn 1993) > 
if the necessary kinematic observations can be obtained as is 
possible for the first time in the M31 shells. Also, there exists 
a simple relation between the surface density of a shell and 
the rate at which stellar mass flows through the shell; when 
appied to the W shelf, this can be used to deduce the orbital 
parameters and the mass of the satellite that produced the 
observed debris. 

Our purpose in this paper is threefold: first, to show 
that this paradigm of a shell debris pattern is consistent 
with the current observations of M31; second, to explain the 
basic physics that gives rise to the spatial and kinematic 
pattern of the debris; and third, to find observational di- 
agnostics to further constrain the model of the debris from 
present and future datasets. In § 2, we give an overview of 
the complex observed structure around M31, and describe 
our scenario to explain a major part of this structure. In 
§ 3, we explain our methods for setting up and running a 
simulation of this scenario. We then compare the simulation 
results to the observed properties of the debris around M31, 
to see whether our model is plausible. In § 4, we discuss the 
physics that governs this model of the satellite debris. We 
consider several observables that can be used to confirm or 
rule out our scenario, and to constrain the physical prop- 
erties of the stream and its progenitor. In § 5, we discuss 
some remaining issues regarding our model, and state our 
conclusions. 



2 SCENARIO OF SHELF FORMATION 

Figure^ shows the density of RGB stars around M31 (from 
llrwin et al .120051) . The sky coordinates ^ and r) in the eastern 
and northern (E and N) directions respectively are "stan- 
dard coordinates" centered on M31. Several features that 
will be important in the discussion below are marked, in- 
cluding the s outhern (S) stream. The overdensity known as 
the NE shelf jFerguson et al.ll2002l: iFerguson etallbOOSD is 
also readily apparent. 

Less obviously, there is a similar enhancement of lower 
surface brightness on the W side, which is apparent on 
close visual inspection in this map and several of its ear- 
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this structure, we display in Figure^ a map obtained from 
Figure la using the Sobel edge-detection operator, which 
turns sharp edges into bright ridges. We fit the NE and W 
shelf ridges in this map with smooth curves, whose interpo- 
lating points are given in Table The fainter W shelf has 
not previously been been singled out as an interesting fea- 
ture in the literature, but it appears clearly enough here to 
investigate its possible origins. 

Several other interesting features show up in Figure 03, 
among them the "Northern Spur", a bright loop at ^ ~ 0.8°, 
rj ~ 1.8°; the "Gl clump" and related structure near the SW 
major axis at the outer regions of the disk; and a shelf in the 
southern stream at r; ~ — 1.0°. Whether these other features 
have a ny relation to the stream itself is not currently known. 

In lPaoer III we found that the southern stream is well 
explained by a baryon-dominated dwarf satellite of mass 
~ 10^ Mq, experiencing its first pass through M31 pericen- 
ter. Generically, in such a collision, much of the satellite's 
mass is stripped and turns into leading and trailing tidal 
streams. These two streams are roughly symmetric about 
the progenitor, though the leading stream spirals in towards 
the primary galaxy and the trailing stream outwards from 
it, as one traces the streams away from the progenitor. In 
M31, the trailing, higher-energy stream is the one that must 
correspond to the observed southern stream, since the sur- 
face density is observed to fall off towards the tip of this 
stream opposite to its direction of motion. The locations of 
the progenitor itself and of the forward, lower-energy stream 
are currently unknown. If the satellite should still be in- 
tact at its next pericentric passage, it again produces two 
streams, which nearly overlay the portions of the original 
stream nearest the satellite. A clear example of such a mul- 
tiple s tream system is shown in lLaw. Johnston, fc Maiewski 
l|2005h using a simulation of the Sagittarius galaxy. 

Since the progenitor is not observed along M3rs south- 
ern stream, which spans a range of ~ 0.4 in orbital phase 
(measured in units of the radial period), it is likely that the 
stream continues at least an equal amount of 0.4 periods 
ahead, taking it past the next pericenter (the second in our 
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Figure 1. Panel a: Map of RGB count density, from llrwin e t al. 12005). The edges of the NE and W shelves that are the focus of this 
paper are marked with red dotted lines. The H13d field of lReitzel et al.. I2006f) and four fields from llbata et al] 120051) discussed below 
are marked with crosses. Panel b: Sobel-filtered version of Panel a, which detects sharp edges in the count map. To create this map, we 
fill in sharp features such as plot annotations and noise spikes with the surrounding pixel values, smooth the map, and apply the Sobel 
operator. We then fit a smooth curve to the shelf edges, which again are shown by dotted lines. 



numbering scheme) and a t least a si Rnificant part of the way 
to the next apocenter. In iPaper ij. we could not trace the 
continuation of the southern stream precisely past the sec- 
ond pericenter where it merges into M31's bulge and disk. 
However, we found that it was constrained to lie on the NE 
side of M31. We found an acceptable (if not optimal) or- 
bital fit putting the progenitor in the area of the NE shelf. 
In contrast, the progenitor would need one more pericentric 
passage (its third) to reach the W shelf, as would any ma- 
terial in the tidal stream. We also see in Figure 1 that the 
surface density of the W shelf is lower than for the NE shelf, 
as if the density were tailing off in the stream forward of the 
progenitor. 

These facts lead us to our basic scenario: as we follow 
the continuation of the southern stream further forward in 
the orbital trajectory, the debris first produces a diffuse shelf 
on the E side of M31, and subsequently a similar shelf on the 
W side. Depending on the current position and state of the 
progenitor, the NE shelf may be made up solely of material 
from the forward continuation of the S stream, or include 
the remnant of the progenitor and the streams created at 
its second pericentric passage as well. In the next section, 
we test this scenario. 



3 SIMULATION OF DEBRIS PATTERN 
VERSUS OBSERVATIONS 

3.1 Method 

In Papers I and II, we discussed the physical geometry of 
M31 and its surrounding system of satellites and streams. 



We then derived a bulge-disk-halo model for M31's grav- 
itational potential (Paper I), and developed a method for 
fitting the orbit of the stream's progenitor given the ob- 
served properties of the stream (Paper II). It is worth repeat- 
ing here M31's heliocentric radial velocity of — 300kms~^ 
\de Vaucouleurs et alJIlQQll) . A lso, we assume the distanc e 
to M31 is Zg = 784 ± 24kpc JStanek fc GarnavichI [l99il . 
which agrees very well with the distance of z„ = 785 ±25 k 



derived from the RGB maps bv lMcConnachie et alJ i2005 
and implies that 1° in angle corresponds to 13.5 kpc. 

To construct a model of the progenitor, we first estim ate 
its orbit, using the methods similar to those in iPaoer III For 
the NE shelf we assume a central position on the sky of ^ = 
1.8°, = 0.7°. We fit a progenitor trajectory to pass in the 
direction (not necessarily the exact location) of this point, 
while still reproducing the observational properties of the 
southern stream, obtained by analytically "stretching" the 
trajectory t o approxi mate the phase-space location of the 
stream (fsee lPaper I J. fo r details)). We also incorporate the 
radii and velocities of the lMerrett et aL (,20C)§) northeastern 
"Stream" PNe into the fit, constraining the orbit to pass 
near their projected radii and line-of-sight velocities. 

The M31 potential model her e is a member of the fam- 
ily of models described in IPaoer"^ and includes a Hernquist 
(1990) bulge, an exponential disk, and a Navarro, Frenk, & 
White (1996, henceforth NFW) halo. The specific parame- 
ters used here imply a low-mass disk. The results are not 
strongly sensitive to the disk mass, though, since the rota- 
tion curve places strong constraints on the overall potential. 

Our approximate analytic method describes only the 
central path of the debris in phase space. To gain a more 
complete picture of the debris pattern, we turn to TV-body 
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Table 2. Simulation Parameters 



Parameter 


Symbol 


Value 


Comments 


Satellite structure 


Progenitor Mass 
Progenitor Radius 
Central Surface Density 
Central Velocity Dispersion 
Particle Mass 
Particle Number 


M, 
as 

O"0,s 

Np 


2.2 X 10^ Mq 

1.03 kpc 

6.57 X 10** Mq kpc-2 
39.1kms-i 
1.68 X 10* Mq 
131072 


Plummer Sphere 
1-dimonsional 


Satellite orbit 


Initial position* 
Initial velocity* 
Final time 


yo 

zo 

VyO 

to 


-34.75 kpc 
19.37 kpc 
-13.99 kpc 
67.34 km s-l 
-26.12kms-i 
13.50kms-l 
840 Myr 




M31 structure 


Total Bulge Mass 
Total Disk Mass 
Total Mass inside 125 kpc 
Virial Mass 

Bulge Scale Radius 

Disk Scale Radius 

Disk Scale Height 

Halo Scale Radius 

NFW profile Virial Radius 


Mb 

M{< 125 kpc) 

^2*00 

M* 
■"^100 

Rb 

Rd 
Zd 

rh 
pt 

^200 
"■100 


3.24 X lO^f Mq 
3.66 X 10^° Mq 
7.3 X 10" Mq 
8.8 X 10" Mq 
9.5 X 10" Mq 
0.61 kpc 
5.40 kpc 
0.60 kpc 
7.63 kpc 
195 kpc 
253 kpc 


Hernquist Sphere 
Exponential Disk 


Disk Central Surface Density 

Halo Density Parameter 

Halo Concentration Parameter 


So 

Sc 

C200 = R20oh'h 


2.0 X 10* Mq kpc-2 

4.41 X 10^ 

25.5 


characteristic density relative to critical 


Maximum Rotation Velocity 




259kms-l 





* X, y, and z here are coordinates in a system where x points to the E and y to the N in the plane of the sky, and 2 points into the sky. 
^ Ma is the mass enclosed within radius R/\ such that the mean density inside is Apc, where A = 100 or 200 (see ^age^T). 



simulations. We continue to use a fixed potential to repre- 
sent the influence of M31, but the progenitor is now rep- 
resented by self-gravitating particles. We consider a purely 
stellar system. It is hard to put much dark matter in the pro- 
genitor without either underproducing the amount of stel- 
lar mass seen in the stream, o r overpro ducing its observed 
width and velocity dispersion jPaoer iJl . The progenitor is 
given a spherical Plummer-law density profile, with rough 
values of the initial mass and scale radius Ms ~ 10^ Mq 
and as ~ Ikpc (later to be refined). We initialize the par- 
ticles of this hot spherical system with the ZENO pack- 
age of J. Barnes. The particle velocities are assumed to 
be isotropic, and are set by solving an Abel integral for 
the energy distribution rather than assuming a Maxwellian 
fsee lBinnev fc Tremainelll987^ . giving an initial configura- 
tion very close to equilibrium. 

Using the orbit from the fit above, the satellite is set in 
motion inbound and slightly past apocenter, minimizing ini- 
tial transients from M31's tidal forces. The formation of the 
progenitor and its accretion into M31's halo are not modeled 
here. Once we determine the current state of the debris, we 



should be able to shed some light on the progenitor's prior 
history, but for now it lies outside the scope of our work. 

We run the simulations with the multistepping tree code 
PKDGRAV tStadel 2 001'). We set the spline softening length 
to a small fraction of the satellite's Plummer scale radius, 
e — O.lfls ~ 100 pc. The gravity tree uses an opening angle 
6 — 0.8, the node forces are expanded to hexadecapole order, 
and the individual particle timesteps Ati are related to the 
particle acceleration at by the criterion Ati < ri{e/ai)^^^ 
with rj — 0.2. 

We ignore the effect of dynamical friction in both our or- 
bital fitting procedure and in the A^'-body simulations. This 
is reasonable given that we are not aiming at an exact fit of 
the observations, since the mass of the progenitor satellite is 
small relative to the total mass of M31 (~10~^), and smce 
the Coulomb logarithm is small as well (In A ~ 3). Even so, 
the specific orbital energy of the progenitor decreases some- 
what during the run, since the energy required to tidally 
disrupt and heat the satellite comes at the expense of the 
satellite's orbital energy. We ignore any perturbations from 
M32 or NGC 205 as well. Their mass is low enough that a 
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significant deflection would require an unlikely close passage, 
and their orbital parameters are not known well enough for 
us to treat this possibility with any reliability. 

To relate the numerical simulations to observable prop- 
erties of the stream, we assume a flxed M/L ratio through- 
out the stellar population, or more precisely a flxed ratio of 
stellar mass to RGB star numbers, since it is the latter class 
of objects that are used to trace the observed stream and re- 
lated debris. We also assume a fl xed ratio of stellar mass to 
PNe number. iMerrett et al.l ll200tf) showed that there is good 
agreement between the number counts of PNe and M31's r- 
band light profile, which traces the old stars and thus the 
stellar mass quite well. 

Our initial orbital fit yielded promising results. We then 
ran several more simulations to choose the mass and size of 
the satellite, essentially through trial and error, to match 
the length of the southern stream and other observable fea- 
tures. We also attempted to optimize the orbital parameters 
further using A^'-body simulations, but without clear further 
improvement. The parameters of our final simulation are 
given in Table |5| We note that we have certainly not ex- 
plored our large parameter space carefully in this search. 
Also, it is not yet possible to make quantitative compari- 
son to some of the most relevant observational data such 
as the RGB spatial distribution, due to thorny issues such 
as incompleteness and contamination from both foreground 
Milky Way dwarf stars and background galaxies. Hence, we 
make no claims that our present satellite orbital parameters 
are optimal. They should, however, suffice to demonstrate 
the basic physical and observational properties we can plau- 
sibly expect from this general scenario. 

Our spherical Plummer model for the satellite with 
initial mass Ms — 2.2 x 10^ Mq and scale radius as = 
1.03 kpc produces a central velocity dispersion of 39 kms~^, 
consistent wit h typical velocity dispersions for galaxies 
of this mass JPekel fc Wod l2003ll . The observed stream 
stars appe ar to have a modera t ely h ig h metallicity, [Fe/H1 



12003 



l2006l: iKalirai 


et al.' '2006tf) 


, and intermediate to old ages 


of 4-11 Gyr (j 


Brown et al. 


,WQQa,,h). Again comparing to 


iDekel & Wod 


ll2003i) ■ this observed metallicity is in very 



good agreement with the assumed progenitor mass. If we 
assume a l/-band M/L ratio of 5 in the progenitor, appro- 
priate for a baryon-dominated galaxy with an old stellar 
population and a Salpeter initial mass function, the central 
surface brightness is Eo « 1.3 x 10* Lq y kpc~^. This is a 
factor of 3 below the mean regression line of Dekel & Woo, 
but is still consistent with the scatter about this line. The 
Plummer law is not particularly concentrated towards the 
center. If we had adopted an R^^''' law for the satellite, the 
central surface brightness would have been about 300 times 
larger, for the same total luminosity and projected half-mass 
radius. 

We run the sateUite model with PKDGRAV for 2 Gyr. 
The first pericentric passage occurs at 0.17 Gyr, giving rise 
to the southern stream. The first apocenter occurs at 0.44 
Gyr and the second pericenter at 0.70 Gyr. We select the 
output that most closely resembles the spatial and kinematic 
properties of the stream and shelf debris. This output occurs 
at 0.84 Gyr past the start of the simulation (Table IJJ. At 
this point the progenitor has been completely disrupted and 
its central stars lie in the location of the observed NE shelf. 



We first tried runs with a small number of particles before 
proceeding to the final run; we did not find any significant 
differences resulting from particle resolution. 



In several of the plots below, we will use only the sim- 
ulated satellite particles for clarity of presentation. How- 
ever, when comparing to the observations, it is essential 
to have an idea of the characteristics of the main stel- 
lar body of M31, which in project ion is coincident with 
muc h of the sate l lite de bris. Both iFerguson et alJ ll2005l) 
and iBrown et all |2006a) found that at 10-20 kpc from 
M31's center, the stellar populations of the spheroid and 
of the stream were quite similar. On the other hand, an 
outer halo component has r ecently been discovered in M31 
iGxihathakurta et alJ 1200,4 llrwin et aLlEoOflt iKahrai et all 
J: IChapman et alJ l200m: incorporating these new ob- 
servations, the mean metallicity of the smooth M31 spheroid 
is seen to decreases radially outward down to [Fe/H] « — 1.4 
at 7? > 60 kpc. In contrast, the stream's mean metallicity is 
more or less constant along its length, out to similar pro- 
jected radii. So depending on location, the characteristics of 
the stellar population may or may not be useful for separat- 
ing the contributions of the stream debris and M31; in some 
regions kinematics may be the only guide. 



Therefore, we construct an Ai'-body particle realization 
of M31 itself to mimic the observed stellar components of 
M31, and combine it with the satellite debris to predict the 
full spatial and velocity pattern. We only add this particle 
realization at the end, rather than e volving it dynamically 
with the A'^-body code. The model in IPaoer ll predated the 
discovery of the metal-poor halo component just mentioned; 
it is important to include this in the model to estimate the 
M31 contribution at large radii. We model the halo stars 
by assuming that a small fraction of the NFW halo mass 
is composed of stars rather than dark matter, equivalent 
to using an NFW profile normalization (5c,* = 350. To fit 
the profile better, we also truncate the Hernquist bulge past 
rtr,b ~ 10 kpc since it would otherwise overproduce the sur- 
face brightness near 20 kpc. We use the default truncation 
scheme in ZENO, which cuts off dM{r)/dr as an exponen- 
tial past the truncation radius. (We neglect this truncation 
in the orbit calculation, since it changes the mass within 125 
kpc by only 0.2%.) With these modifications, the composite 
stellar profile exhibits a change in slope at ~ 30 kpc similar 
to that in the observed profile iGuhathakurta et al.. .20051 : 
llrwin et aLlEjoilKaliran et alJl2006ah. 



For our TV-body realization of M31, we use the same 
particle mass as for the progenitor particles, and again 
use the same constant scaling of RGB or PNe number to 
stellar mass. Given that there are slight but observable 
stellar populati on differences betwee n M31's giant stream 
and spheroid l| F^rgusqn ^l.| [2(^(^)51: [ Brown et alJ l2006al : 
iKalirai et alj lioOGbjatlBroTOL^anlioOG^ constant 
scaling assumption is a slight oversimplification. Moreover, 
there are significant observed deviations from our simple 
model in the form of spiral arms, warping of the disk, and 
elongation of the bulge and inner disk. Our model should 
thus be taken as an illustrative proxy for the M31 stellar 
component, rather than a precise description. 
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Figure 2. The spatial distribution of satellite debris particles in the simulation, centered on M31. Panel a: sky projection in standard 
coordinates. This can be compared to the observed star-count map in Figure la. The red dotted lines show the observed shelf edges from 
that figure. For clarity, a randomly chosen subset of particles is shown in this and most subsequent figures. The particles of our static 
M31 model are also omitted; the disk component of that model would dominate much of the NE shelf but only a small portion of the W 
shelf. Panel b: the view of the debris looking down the north axis in Panel (a). "Depth" dM3l increases with distance from Earth. The 
scale on that axis matches the angular scale at M31's mean distance of f» 780 kpc. The NE shelf is directed mainly toward us, while the 
W shelf is oriented away from us close to the plane of the sky. 



3.2 Shelf morphology 

We begin the comparison to observations by examining the 
pattern on the sky of the NE and W shelves. The spatial 
distribution of the particles in our simulation is shown in 
Figure 121 The debris resembles the observed shelves (shown 
in Figure 1) on both the E and W sides, with a much larger 
surface density on the E side. The edges of these shelves 
are at similar, but in some places slightly smaller radii to 
the observed edges. The shelves also cover about the right 
azimuthal range, though the faintness of the features and in- 
terference from M31 itself make it difficult to make a precise 
comparison 

We note that the projected radius of the NE shelf from 
M31's center, as shown in Figure and given in Tabled 
appears to vary with azimuth, instead of forming a circu- 
lar arc. This may be a significant clue to the nature of the 
shelf. In contrast, the edge of the W shelf appears more cir- 
cular. The simulations appear to reproduce the behavior of 
projected radius as a function of azimuth for both shelves. 

We next turn to the s urface densities. In the star-count 
map of llrwin et al.l (|200^, the brightness of the NE shelf 
appears similar to the denser, innermost distinct parts of 
the S stream (at Rproj ~ 1-6°), while the brightness of the 
W shelf appears more similar to the more tenuous middle 
parts of the S stream (at Rproj ~ 2.5°). To obtain a crude 
estimate of the surface densities in the latter regions, we first 
construct a simple estimate for the southern stream surface 
brightness using the trans verse and longitudinal profiles in 
iMcConnachie et alJ i2003h . We normalize these profiles by 
using the stream's V magnitude of My = —14.0 or V lumi- 
nosity of 3.4 X L r^ ^, in the range 15 kpc < Rproj ^ 40 kpc 
from M31's center lllbata et alJl200]]) . Assuming a stellar 



M/Lv ~ 5 as before, this gives a stream stellar mass within 
this radial range of ~ 2 x 10* Mq. Computing the surface 
density of the resulting mass model in the areas mentioned 
above, we can then estimate that the NE and W shelf sur- 
face densities are very roughly ~ 1.0 x 10® Mq kpc~^ and 
~ 8 X 10^ Mq kpc"'^ respectively. We arbitrarily define the 
inner radial boundary of both shelves to be 60% of their 
outer radii (listed in Table in order to isolate the regions 
in which the shelves have the highest contrast against the 
smooth M31 components; the NE shelf has an area 150 kpc^ 
and the W shelf has an area 200 kpc^. This then implies 
masses of ~ 1.5 x 10** Mq within both shelf boundaries. We 
have not attempted here to subtract the background density 
of M31 disk and spheroid stars, so these surface densities and 
masses are only lower limits. Also, these estimates only ap- 
ply to the mass inside our assigned boundary; if the shelf 
mass distribution continues radially inwards onto the face 
of M31, the total mass could be much larger. Given all the 
uncertainties, these are clearly only order-of-magnitude es- 
timates, but even so should be good enough to check the 
plausibility of the model. 

We compute the mass of the satellite debris particles in 
the simulation that lie within the boundaries of the NE and 
W shelf regions as defined above, and find 1.2 x 10* Mq 
for the NE shelf and 1.4 x 10* Mq for the W shelf (al- 
though these regions clearly do not select all the mass in 
the NE and W radial lobes). These imply surface stellar 
mass densities of 8 x 10^ Mq kpc~^ for the NE shelf and 
7 X 10^ Mq kpc~^ for the W shelf. The former surface den- 
sity in particular is a slight underestimate because it does 
not extend all the way to the observed shelf edge. Inclusion 
of the particles in the static M31 model raises these sur- 
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face densities to 1.1 x 10® Mq kpc^^ for the NE shelf and 
1.0 X lO^MQkpc"^ for the W shelf. Comparison of these 
values with the rough observational estimates above shows 
that our model has passed another test. 

M31 is close enough that the di stribution i n thre e spa- 
tial dimensions is measurable. Ferguson et al.l ll2005ll used 
color-magnitude diagrams with HST/ACS to measure the 
brightness of red clump stars in the stream, 21 kpc in pro- 
jection from the center of M31, and in the NE shelf, at 
£, = 1.26°, 77 = 0.24°. They found that stars in the "NE 
shelf" fields indeed lay closer to us than those in the "Giant 
Stream" fields, by a factor 1.07 ± 0.01. With the "depth" 
d (defined here as the distance behind M31) in the stream 
field es timated to be (50 ± 20) kp c by interpolating the re- 
sults of lMcConnachie et alJ (l200,'^) . thi s would put the shel f 
at depth dpos = (0±20) kpc. However. ISrown et al.l (l2006al) 
measured the red clump brightnesses at a nearly identical 
point in the stream and in a minor-axis, spheroid-dominated 
field, and found the stream there was only (11 ±5) kpc fur- 
ther than the spheroid position, which is likely at d ~ 0.^ 
From this estimate we can infer the shelf is at d « —40 kpc, 
in front of M31. 

Our model implies the material forming the NE shelf 
hes on average well in front of the center of M31, and hence 
even farther i n front of the southern stream (see Figure I^J. 
In the field of iFerguson et all (|200^, the median depth of 
the simulation particles is d pos = —31 kpc, in a greement 
with the estimate based on iBrown et alT (12006311 . though 
the distribution is highly skewed with a small second peak 
at « —5 kpc. The model thus appears to pass the test of the 
observed distances as well. 

3.3 Velocity measurements 

The overall velocity pattern in our simulation is shown in 
Figure |3] We will discuss the physical origin of the features 
in this plot later, in § 4.3. Here we only wish to compare the 
simulated velocities to the observed ones to see if they agree 
in their main respects. At present, there are three large- 
scale surveys that can constrain the overall velocity pat- 
tern of the ob s erved debris: the wide-area survey of PNe by 
iMerrett et al.l i2003l 2006) and two ongoi ng surveys of RGB 
velocities in selected widely spaced fiel ds (|lbMaL£Laljl20Q4|; 
Ibata et al. 2005; Chaoman ct al. 2006; Guhat hakurta et alJ 
2006i:.Kah rai et al. 2006b a; Gilb ert et al 200^). 

The PNe detected by Merrett et al. include a component 
that moves opposite to the disk rotation on the NE side, 
with a narrow velocity width, suggestive of a continuation 
of the southern stream. The PNe indicated by Merrett et 
al. as "Stream" or "Stream?" are shown as triangles in the 
right panel of Figure |3] This apparent stream is unlikely 
to be produced by the usual stellar components of M31, so 
it is strongly indicative of some kind of tidal debris. The 
simulation produces a similar "stream" of material, due to 

^ Flattening of the spheroid could introduce a systematic offset 
to this estimate, to larger distances if the spheroid is aligned with 
M31's disk. Howe ver, the consistency to within < 0.03 mag of the 
spheroid fields of [Bro wn et alJ J20063), at Rproj = 12 kpc, and 
iFerguson erall <2005h . further out at Rproj = 20 kpc, suggests 
that the Brown et al. field is shifted by < 0.03 mag or < 12 kpc 
with respect to M31's nucleus. 



a set of stars in the N part of the NE shelf which have similar 
velocities; like the observed PNe stream, these overlap the 
major axis of M31, as shown in Figure |3| below. 

The Merrett et al. survey extends into an irregular area 
in the halo region of M31. In this region, Merrett et al. used 
wide-field images to optimize the choice of fields for closer 
study. These fields are sparse, so the survey coverage may 
change with position; we will ignore this potential compli- 
cation. Figure 0^ shows the survey results in the area of the 
W shelf; the survey coverage within the W shelf boundary 
(dashed red line) and between the dashed blue radial lines 
is almost complete, while it is sparser but not zero outside 
the W shelf boundary. We exclude the PNe that Merrett et 
al. find to be associated with NGC 205 on the basis of their 
positions and velocities. The same plot shows the simulation 
particles, including both the satellite debris (in green) and 
the particles from the static M31 model (in black). Although 
the number of PNe is limited, it seems that the density of 
PNe does not fall off all the way towards the W shelf; be- 
yond R ~ 0.7° it remains roughly constant, suggesting that 
the PNe also trace the shelf. We plot the velocities of the 
PNe and simulation particles within the dashed selection 
boundary in FigureUJ). The satellite debris in the simulation 
forms a triangular shape in this plot. Interestingly, the PNe 
at large Rproj mostly fall near the boundary of the simulated 
satellite debris; it appears that increasing the shelf radius or 
particle velocities slightly in the simulation would result in 
even better agreement. 

To test this apparent agreement quantitatively, we com- 
pare two hypotheses: that the distribution matches the pure 
M31 particle distribution, or that it matches the sum of the 
M31 and satellite particles. We bin the particles in a grid 
covering the relevant area in projected radius and velocity, 
after convolving their velocities with a Gaussian distribution 
of 17kms~^, which is the random error in the PNe veloc- 
ities. We then normalize the distributions, both including 
and excluding the satellite debris particles, to unity. The 
product of these binned distributions at the locations of the 
PNe then gives the likelihood function of these two hypothe- 
ses. We find that the likelihood ratio is about 180, meaning 
that the model with satellite debris is 180 times more likely 
than the model that excludes it, assuming equal prior prob- 
abilities. Furthermore, if as suggested by Figure 21 we add 
30kms~^ to the velocities of the satellite particles as might 
be induced by a change to the overall angular momentum, 
the likelihood ratio increases to 4500. Thus these tests ap- 
pear to show that some of these PNe are bona-fide members 
of the W shelf, and confirm its overall kinematic structure. 

The remaining M31 radial velocity surveys use RGB 
stars. Previous analyses of substructure in these surveys 
mostly focus on the giant southern str e am, the extended 
disk-like com ponent Jibata et al.l l2004t llbata et al] l2005l : 
(Guhathakurta et alj| 2006|). and"tlie loop- like features near 
NGC 205 lIMcConnachie et aljliooi) . but there are some 
hints of other structures. The stream of "counterrotating" 
PNe may have a counterpart i n the inner disk RGB velocity 
measurements of llbata et al.l ll2005l their Figure 9) , though 
this fe ature is not discussed in that paper. Chapma n et alJ 
i2006li mention kinematic substructure in the NE that they 
attribute to a continuation of the southern stream, but the 
exact location of the responsible fields is not clear. 

The RGB surveys obtain a much higher density of tar- 
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Major-axis position X^^j, (degrees) Major-axis position X„3j (degrees) 



Figure 3. Plots of the velocities of satellite debris particles relative to M31 versus major axis position X, split on the basis of minor-axis 
position Y into two panels. X increases toward the SW, and Y toward the NW. Left panel: southeast portion {Y < —0.3°). The dense 
concentration of po ints at the lower right edge is the southern stream. Right panel: central and northwest portion {Y > —0.3°). The PNe 
in M31 notated bv lMerrett et aP i200(i) as "Stream" or "Stream?" are marked by red diamonds. The dense concentration of Af-body 
particles within 60kms~^ of the lower left edge are marked in green; we regard these particles as an analogue to the Merrett et al. "PNe 
stream" . 
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Figure 4. Plots of Af-body particles and observed PNe in the W shelf region. Objects with minor axis coordinate Y < +0.25° are 
omitted. Panel a: sky position in M31-axis coordinates X and Y. Green points s how satellite debris p articles, while black points show 
particles from the model of M31 itself. Orange triangles show observed PNe from lMerrett et al.l (jOQCj), excluding those noted there as 
associated with NGC 205. The red dashed line shows the observed W shelf boundary from Figure^and Table^ The sky coverage in the 
PNe survey extends roughly out to the outer blue dashed arc. Only objects shown within the blue dashed lines are selected for the next 
panel. Panel b: projected radius versus line-of-sight velocity. The symbols are the same as before. The PNe velocities have a measurement 
error of ITkms"'^. Note the group of PNe close to the positions and velocities of the simulated shelf particles; this association is tested 
statistically in the text. 
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Figure 5. Th e velocity distri bution relative to M31 in the 
field H I 3d of Rcitzol ot al. 1 2006) (also discussed in*K alirai et alj 
I2OO6IJ) . The field used to select the simulation particles is taken 
to be a square 0.3° on a side, to ensure a large sample of particles. 
Simulation results using the satellite debris alone are shown by 
a red dashed line; results using the M31 galaxy model alone are 
shown by a green dotted line; and results using both satellite and 
M31 particles are shown by a thick solid line. The distribution 
of observed RGB star velocities is shown with a cyan histogram; 
this is normalized to have the same total density as the full sim- 
ulation histogram . The disk velocity from the empirical model in 
llbata et alj l2005l) is shown by an arrow. The observations show 
at least a qualitative match to both the prograde (disk) and retro- 
grade (satellite debris) components in the simulation, as discussed 
further in the text. 



get objects than the PNe survey, albeit with a much sparser 
spatial coverage. Thus it is worth plotting the velocity distri- 
bution in each fiel d as a separate histogr am. We begin with 
the fi eld H13d of iKahrai et all 1120061^) and iReitzel et all 
i2006i) . The observational and simulation results are shown 
in Figure We include three curves that show the re- 
sults of the simulation for the satellite debris alone, for the 
M31 model stars alone, and for the combination of the two. 
Here we assume the total mass density in the observations 
matches the total found in our simulation, and normalize the 
observational histogram accordingly. This method allows us 
to avoid the complex issues of the selection and detection 
efficiencies in the spectroscopic survey. 

A substantial population at negative velocities, oppo- 
site to the sense of disk rotation, is evident in the observa- 
tional histogram. This component is not visible at all in the 
galaxy model, but it matches the velocity distribution of the 
satellite debris quite well, even if the overall normalization is 
somewhat off. Another interesting feature is the difference in 
mean velocity, velocity dispersion, and possibly total mass 
between the observed and predicted disk components; we 
will return to this point below. 

The spectroscopic survey of llbata et alJ (|200^ also in- 
cludes 16 fields in the outer disk. In at least 7 of these fields, 
the velocity distribution shows a strong, narrow component 
near the expected velocity of the disk, which forms the focus 



of that paper. None of the other fields is dominated by the 
disk component, although a disk component emerges when 
the fields are stacked together. As shown in Figure four 
of these latter fields lie in the regions possibly dominated by 
debris associated with the NE and W shelves, namely FIO, 
Fll, and F13 in the NE shelf, and F4 in the W shelf. The ob- 
served velocity distributions in these fields are shown in Fig- 
ure|ni together with the prediction from our simulation. One 
complicating factor is that Ibata et al. remove all stars with 
velocities greater than 200kms~^ (or — lOOkms"^ heliocen- 
tric), to reduce contamination from MW dwarf stars. They 
also plot histograms of the velocity relative to their disk ve- 
locity model; we have transformed these to M31-centric ve- 
locity using their model velocity in the center of their field, 
but this correction neglects the slight variation of this model 
velocity across the field. The shape of the histogram is thus 
not exact, especially near the velocity limit of 200 km s^^ 
(shown with a dotted line) 

In F4, FIO, Fll, and F13, stars are detected in a broad 
velocity distribution up to |«| ~ 200kms~^ relative to M31, 
as shown by the light cyan histograms in Figure |S| A glance 
at Figure |S] shows that in fields FIO and Fll, the predicted 
satellite and static M31 velocity distributions are actually 
fairly similar. Thus, the velocity distribution in these two 
fields cannot serve as evidence for or against our scenario 
without a deeper understanding of the slight differences. 
The histograms from the simulations suggest that the satel- 
lite debris dominates the distribution in both fields; dif- 
ferences in the stellar populations might serve to test this 
result, though this also could be dif ficult as these differ- 
ences are quite subtle at this radius jFereuson et al.ll2005l; 
Brown et al. 2006a b?). 

In field F4, the strongly bimodal satellite distribution 
biased to positive velocities improves the agreement between 
the simulated and observed velocities, though the observed 
distribution extends to even higher velocities than shown by 
the simulated satellite debris. It is worth noting that the 
problem of foreground MW dwarf contamination increases 
with the velocity. The observed and simulated results in field 
F13 are similar to those in H13d above; again there is a 
population of observed stars at negative velocities, which 
are explained fairly well with the satellite debris. These two 
fields are thus encouraging for the debris model. 

Field F13 also resembles Field H13d in that there is a 
strong predicted galaxy disk peak that is not seen in the 
observations. Once again there may be a smaller disk peak 
shifted to lower velocity by 40-60 km s~^. The velocity cut 
at 200kms~^ imposed by Ibata et al. to select M31 stars 
could be responsible for at least part of the difference, as 
it would remove many (though not all) disk stars in our 
Ai'-body model for M31. Moreover, the warp of MSl's disk, 
which is quite strong this far along the NE major axis, may 
cause the disk contribution to be substantially lower in field 
F13 than in H13d (F13 being slightly off the major axis 
in the direction away from the warp). However, the H I in 
this r egion shows a peak at 220kms~^ iNewton &: Emerson! 
(^73), exactly as predicted by the disk model, even though 
the H I disk itself is warped. A further possibility is that 
unlike the H I, the stellar velocities are perturbed from cir - 
cular orbits. The extended disk found bv llbata et alJ (|200^ 
shows large velocity and density fluctuations relative to a 
smooth disk, and the results in fields H13d and F13 may be 
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Figure 6. Same as Figure|^for the four spectroscopic fields observed for RGB star velocities bv llbata et alJ (2223) that lie in the shelf 
regions. Otherwise the caption to the prev ious figure applies here. The observational histograms are generated by taking the histogram 
of "disk lag velocity" in llbata et al.l j20o3) and shifting by the velocity in their model of the disk at the center of the field. This neglects 
the slight variation of the disk velocity across the field. In these observational histograms, stars with positive velocities greater than 
200kms-l relative to M31 are removed to reduce Galactic contamination, as indicated by the dotted line. The observations show evidence 
for satellite debris in at least Panels (a) and (d), as discussed in the text. 



a further example of these fluctuations. A difference in H I 
and stellar velocities might be expected if the extended stel- 
lar disk resulted from an accreting satellite, as in the model 
of (Pena rrubia 2006). 

In summary, it appears that the debris model studied 
here explains a number of irregularities seen in the spatial 
and velocity distribution of the stars and PNe in and around 
M31. The strongest pieces of evidence for the model are the 
shape of the NE shelf on the sky, the apparent stream of 
counterrotating PNe and RGB stars seen along the NE ma- 
jor axis, and the agreement of the W shelf in the model with 
the spatial pattern seen in RGB stars and the kinematic dis- 
tribution of PNe in the same area. The measured distance 



to the NE shelf stars is also consistent with our model. The 
fields observed for RGB velocities to date are not in re- 
gions that are optimal for testing the model, even though 
the existing fields are encouraging for the debris model, so 
we can anticipate great improvement in the observational 
constraints on this model. We now examine how future ob- 
servations can test the model for the debris pattern, how we 
can understand its origin physically, and how we can use the 
debris to better understand M31 and its environment. 
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4 PHYSICAL INTERPRETATION AND 
DIAGNOSTICS 

4.1 Model for collisionless shell systems 

The NE and W shelves in our model appear similar to the 
"shells" observed around a significant fraction of elliptical 
and SO galaxies. One mechanism for forming s hells involves 
the disruption and folding of a thin, cold disk ('Quinn'19841; 
[HcrriQuist & Quinn 1988; Barnes & Hcrnauist 1992). This 
scenario may well occur in nature, but it would tend to pro- 
duce sharp, asymmetric features with no clear orie ntation 
aroun d the galaxy center (for example, see Figure 3of|n umnl 
Il984l) . and hence be a poor fit to the best-defined examples 
of shell patterns in ellipticals, as well as the tangentially ori- 
ented shelves discussed here. A more likely creation mech- 
anism in our case is t he accretion of a satellite on a nearly 
radia l tra.iectory (e.^.. | Schweizeilll980l: iHernauist fc QuinnI 
Il988l: Ib arnes fc Hernauistlll992H . Once the satellite is dis- 
rupted in whole or in part by its first close passage, the 
correlation between orbital energy and orbital period sorts 
the newly unbound stars in order of energy into a dynami- 
cally cold stream, with the lowest-energy stars executing the 
smallest and fastest orbits. The cooling of the stream means 
the progenitor need not itself be dynamically cold to produce 
a coherent shell, though this is possible as well. Subsequent 
passages near the center of the galaxy, combined with the 
dispersion in the angular momenta of the stars, disperse the 
stream stars in a variety of directions. The energetic sorting 
ensures that they nevertheless reach turnaround at a simi- 
lar radius at any given time. This ongoing phase- wrapping 
process creates a series of caustic or "fold" features located 
on spherical shells, whose radii decrease along the orbit. 

This scenario is much like our curre nt understanding 
of the creation of M31 's southern stream ilbata et al ] |200ll 
l2004l:lFont et al.l2006l : Paper I). Hence shell-like features are 
expected in M31 as long as the debris extends far enough 
ahead of the stream. The previous section provides evidence 
that this is indeed the case. In this subsection we discuss in 
more depth the physical properties of radial shells. 



4-1.1 Stellar orbits in shell systems 

To provide a concrete and analytically tractable example, 
let us make a number of approximations. First, we as- 
sume a power-l aw potent ial for the primary galaxy with 
$(r) oc — r'^. In iPaper III we found this was a reasonable 
approximation to the actual potential of M31 in the radial 
range 3-100 kpc, with k ~ —0.4. 

When the satellite collides with the primary galaxy, the 
tidal forces alter the energies of stars depending on their 
positions and orbits within the satellite, much like the grav- 
itational slingshot of a planetary probe. In turn, the dis- 
persion in their specific energies is related to the progenitor 
satelli te's mass Ms and pericenter rp, as fo und in|johnstonl 
il998h and confimed by the simulations in lPaper ij . Using 
the latter simulations we estimate 



aE ~ 0.7V, (rp) 



GMs 



rpVcf{rp) 



1/3 



(1) 



which we also found to be valid for our run in this paper. 
Here Vc{r) = {GM{r) /rY^^ is the circular velocity. The dis- 
persion then increases with decreasing rp. The combined 



requirements of a large dispersion in the energy and large 
angular deflections at pericenter explains why shell systems 
are most easily modeled using satellites with small rp, i.e., on 
nearly radial orbits. We should note that we have not tested 
the validity of EquationQoutside the regime of nearly radial 
orbits and an M31-like potential. 

Next, we assume the newly liberated stars travel on 
purely radial orbits, i.e., we ignore the angular momentum. 
This is sufficiently accurate in our scenario, where the un- 
perturbed orbit of our progenitor has a pericenter 45 kpc 
and apocenter 2 kpc. The stars then take on orbits that 
are described by the energy alone, or equivalently their or- 
bital periods tr- Let the turnaround radius of the progen- 
itor be rto and the circular velocity there be Vto- Then 
K(r) = Vto{r/rtof'^. 

For an arbitrary star, let us define a time scaling factor r 
proportional to the orbital period of the star, r = tr/tro- The 
self-similar nature of the orbit family leads to the scalings 
rt oc r^/t^-fe), Vt oc r'^/t^"''-), and E oc r^'^/f^-'^). Stars with 
lower energies move faster, as long as fe < 0. The radius and 
velocity of a star with a given r can be obtained from the 
orbital trajectory of the progenitor using the relations 

r{t) = r2/(2-'=Vo(t/T) , (2) 

v{t) = r'=/(2-'='t;o(t/r) . (3) 

Here we are treating the progenitor as a test particle, i.e., 
we ignore any change in the progenitor's orbit from tidal 
disruption and dynamical friction. The exact orbital period 
of the star is 



B[l/2,(2-fc)/(-2fc)] 



Vt 



(4) 



where B is the beta function. 

The turnaround radius of a shell is found for dr/dr — 0, 
which results in the requirement 



t 



voit/r) 



-ro{t/r) = 



(5) 



solving for the roots of this equation yields a series of spher- 
ical shells at fixed orbital phases. The radius of the shell at a 
given orbital phase is then = r^/(^~*Vo(i/T) oc t^/^^"*^). 
The shell moves outwards with time, as stars with higher 
and higher energies arrive at the shell orbital phase. The ra- 
tio of the radii of two different shells is constant with time. 
These results are specific to the power-law potential, but 
also give a sense of the shell properties in the more general 
case. 



4-1.2 Idealized surface density 

To understand the sky pattern of a shell system, we must 
realize that the stars in a given shell may cover only a small 
fraction of the sphere, particularly early in its history as 
is the case for our simulation (cf. Figure EJ . The observed 
properties of the shell depend on whether or not the shell 
surface lies tangent to the line of sight anywhere. These two 
cases are illustrated in Figure 

If the shell surface is tangent to the line of sight, as in 
the leftmost shell in the figure, then the caustic feature at 
turnaround is visible in projection as a circular boundary of 
the shell material. Let us here neglect locally the variation of 
the energy and mass flow rate along the stream. All stars in 
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such a stream share the same shell radius . Near this shell 
radius, the velocity of a star as a function of time t is given 
hydr/dt^-r-^V^{rs){t-t,) = -%/2K(r,)[(r, -r)/r«]i/2, 
where r(ts) = and 14 is the circular velocity. Combining 
the inflowing and outflowing stars, we use this to derive the 
mass density of stars as a function of radius: 



p(r) 



dM 



1 



47rr-2|dr/dt| dt 23/% /n Vc(r,) 



/r^_-r\-i/2 dM 



Here dM/dt is the rate at which mass flows through the 
stream (for example, the amount of mass that reaches apoc- 
enter per unit time), /n is the fraction of the total An sr of 
the sphere that the shell covers; put another way, it specifies 
the ratio of the local mass fiow rate per unit solid angle to 
its average over the sphere. 

When integrated over constant projected radius, this 
gives the result that the mass surface density S approaches 
a constant just inside the shell radius. 



nrs) = 



fe 



dM 



2fnrsVc{rs) dt 



(6) 



and drops to zero outside of it. Here fe represents the frac- 
tion of the mass that is present along the line of sight com- 
pared to a full spherical shell; at the shell radius fg — 1 
(otherwise it would not even be visible), but it can drop 
off at smaller radii, as can be seen in the leftmost shell in 
Figure |7| This factor partially compensates for the solid an- 
gle factor /q. The reason for the constant surface density 
at the edge is that as the line of sight moves inward from 
the shell edge, it reaches a smaller minimum radius, but also 
takes a shorter (increasingly vertical) path through the ma- 
terial at the edge of the shell; these two effects cancel each 
other out. The relationship just derived allows us to obtain 
the local mass flow rate from the surfa ce brightness of the 
shell, which we make use of in § 4.4. jHernauist fc QuinnI 
119881 previously noted that the surface density approaches a 
constant, without relating it to dynamical quantities.) 

If the shell surface is not tangent anywhere to the line of 
sight, as in the rightmost shell in Figure [7| then the extent 
of the shell in projection is limited by the angular coverage 
rather than the shell radius. The dropoff outside the shell 
region need not be sharp, and the boundary of the shell on 
the sky need not be a circular arc. The surface density in 
this case can be estimated from the average over the sphere. 
A single spherical shell corresponds to a single orbit, so the 
total mass contained within the shell equals trdM/dt, where 
tr is given by Equation 3] The total surface area of a full 
spherical shell on the sky is Trr^, so the average surface den- 
sity is 



(-fe) 



B[l/2,(2-fe)/(-2fe)] 



fo 



dM 



in Ts Vc{rs) dt 



(7) 



where we have included corrections for the angular cover- 
age as before. For fc — —0.4, the constant in front is 0.76, 
slightly exceeding the value of 0.5 obtained at a true shell 
edge (c.f. Equation , although typically the value of fg 
will also be smaller than at an edge. Thus there is no strong 
limb brightening or darkening; the surface density does not 
depend very much on whether the line of sight is at a true 
edge or not. 
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Figure 7. Idealized illustration of partial shells in a spherical 
potential. Two shells are shown, which are identical except in their 
orientation. Each shell fills a cone with an opening angle of 30° ; 
the radial distribution assumes a steady stream in a power-law 
potential. The leftmost shell is tangent to the line of sight, and 
thus its outer boundary on the sky is defined by the shell radius. 
Note the dark region where the shell's "end cap" is visible. The 
rightmost shell is not tangent to the line of sight, and as a result 
has a less distinct edge, and a boundary that varies in projected 
radius. 



4.1.3 Idealized velocity pattern 

iMerrifield fc KuiikenI lll998l henceforth MK) studied the 
kinematics of shells. They assumed a shell composed of stars 
on monoenergetic, low-angular-momentum orbits lying tan- 
gent to the line of sight. They then found that these stars 
have a characteristic kinematic signature, which depends 
only on the shell radius and the total mass contained within 
it. In the space of line-of-sight velocity vios versus projected 
radius R, the MK result is that stars near the shell radius 
inhabit a region with a triangular boundary: 

R 



\vioa\ < Vc(r-a) 



(8) 



The stars fill this region unevenly, congregating at the 
boundary. This boundary is again a caustic feature, which 
forms because the velocity component in the line-of-sight 
direction z goes to zero at both z — and r = rs, so it 
must reach an extremum at some intermediate point along 
the line of sight. 

For the discussion below, it is important to understand 
the relationship between the position and velocity of the 
stars. We illustrate this in Figure|Hlfor a full spherical shell. 
The region in physical space that corresponds to the velocity 
caustic is marked by a red line. Caustics occur at positive 
velocities in two sets of stars: outbound stars lying more 
distant than the galaxy center, and inbound stars lying in 
front of the galaxy center. There are negative-velocity caus- 
tics made up of the other two combinations. 

In the case considered by MK, the velocity distribution 
is symmetric around zero, so that the slope is equal and 
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Figure 8. Location of the velocity caustics, in an idealized model 
of a spherical shell of stars. Here we assume a spherical, power- 
law potential with $ oc — r"" *, and a monoenergetic stream. 
Projected radius Rproj and line-of-sight distance d are plotted in 
units of the shell radius. We have reflected the pattern around 
the vertical axis. Since the shell geometry is spherically symmet- 
ric, rotating the entire figure around the vertical axis gives the 
full three-dimensional pattern. The shell boundary is shown by 
the dashed line. The contours show the the line-of-sight velocity 
relative to the maximum velocity along each line of sight, i.e., 
at a fixed Rproj', contour values are at (—1,-0.8,— 0.6. ..0.6, 0.8,1). 
The maximum and minimum velocities along each line of sight 
occur on the caustics marked with red lines, which are labeled 
appropriately for the outgoing stars. For the incoming stars the 
labels would be reversed. 

opposite on the positive and negative velocity sides, and 
vios = at rs. Also, there is a degeneracy between the stars 
on the near and far sides of the galaxy; outbound stars with 
z < have the same sign and distribution in velocity as in- 
bound stars with z > 0, and vice versa. Several effects may 
break this degeneracy in a more realistic case, as discussed 
below. The method suggested by MK has not yet been ap- 
plied to measure the potential in observed shell galaxies, due 
to the low surface brightness of the shells and the difficulty 
in estimating the bimodal velocity signature from weak ab- 
sorption lines. 



4.2 3-D geometry of the debris in M31 

In the discussion above, we have emphasized the role played 
by the angular coverage of the shells. This is important in 
our model since the shell features cover nowhere near the 
full extent of the sphere, as can be seen in Figure |5| It is 
worth taking a closer look at the origin and effects of the 
spatial structure. 

In our simulations, the NE shelf is made up of both 
material from the leading stream from the first pericentric 
passage, which is spread out into a fiattened structure at 
the second pericentric passage, and material stripped from 
the remnant of the progenitor at its second pericentric pas- 



sage, when it completely disrupts. The disruption is largely 
a product of the central density chosen in our initial condi- 
tions, since a more concentrated satellite survived the same 
passage in one of our test runs. While complete disruption is 
in accord with the lack of an obvious visible progenitor, it is 
not required in order to reproduce the shelves. The forward 
stream is further distorted at the third pericentric passage 
into a conical structure, which forms the W shelf. 

The solid angle covered by the simulated NE shelf debris 
is actually even more limited than is apparent in Figure |5| 
When viewed in three dimensions, it resembles a pair of fans 
joined at the edges, corresponding to the outbound and in- 
bound sections; each of these fan layers is only a few kpc 
thick, as shown in the second panel of Figure |5| Because 
this portion of the debris shell does not lie tangent to the 
line of sight, its boundary viewed in projection is defined 
by the angular coverage and not the shell radius, and thus 
varies in projected radius as a function of azimuth. The lat- 
ter property also seems to be shared by the observed NE 
shelf (see Figure and Table , although further work is 
warranted to confirm this result. 

The southern stream is observed out to near its 
turnaround radius, so its length scale is well constrained. 
Above we found that the radii of different shells maintain a 
fixed ratio with time. Since the southern stream is in essence 
the shell preceding the NE shelf, this suggests that the shell 
radii are highly robust for any model that fits the stream. 
Indeed, we found it difhcult to change the shell radii signif- 
icantly for simulations with acceptable stream properties, 
despite a variety of initial conditions and orbital phases. 

The outbound fan on the E side overlaps the region of 
the velocity caustic, as seen in Figure |3] This leads to a new 
interpretation of the Merrett et al. PNe "stream": it does 
not represent an isolated stream of stars traveling along a 
single path, but rather represents an intersection between 
an approximately two-dimensional distribution of stars and 
the mainly two-dimensional region of space where the veloc- 
ity caustic occurs. The various panels of Figure |3] indicate 
in green the "PNe stream" simulation particles previously 
marked in Figure |3 The red symbols in the last panel show 
the location on the sky of the "stream" PNe found by Mer- 
rett et al. , which agree well with the location of our apparent 
stream. From the top view of the structure (Figure|3J;), how- 
ever, it can be seen that the "PNe stream" particles do not 
form an isolated density concentration, and actually inhabit 
a slightly curved region on the fan structure rather than fol- 
lowing a single orbit. The "PNe stream" stars therefore just 
constitute a portion of a larger structure, singled out by the 
kinematics in projection along the line of sight. 

The inbound portion of the fan has its own, inbound 
"stream", which is offset about 0.5° to the SW of the out- 
bound "stream" . This unfortunately mainly lies just outside 
the region surveyed by Merrett et al. for X < —1.0°, and 
for X > —1.0° the velocities mainly lie within the range 
expected from disk stars. However, this returning stream 
should be detectable in a wider PNe survey or a sufficiently 
large RGB sample. A sample near ^ « 1.7°, ~ 0.7° may 
separate the returning component from the disk spatially, 
while a sample near ^ « 1.6°, rj ~ 1.2° may separate it kine- 
matically as it should have a lower velocity there than that 
of the disk. The existence of a returning stream is a critical 
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Figure 9. The NE shelf's fan structure. In each panel, the dots indicate A'^-body particles; those designated as "Merrett PNe stream" 
particles in Figure |2l are marked in green. Panel a: the top view of the structure. M31's center is located on the right at (0,0). The 
horizontal axis points along the main axis of the structure away from M31, while the vertical axis points across the other elongated 
direction of the structure. Panel b: a side view of the NE shelf structure, at right angles to the previous panel. This clearly shows the 
outbound (top) and inbound (bottom) layers of the double fan structure. Again, M31's center is on the right at coordinates (0, 0). Panel 
c: distance d relative to M31, versus projected radius on the plane of the sky. Only those particles in the NE half of M31 {X < 0) are 
shown. Compare the position of the "stream" particles here to the spatial region occupied by the velocity caustic in Figure IHl Panel d: 
th e sky view of the NE shelf structure, plotted in M31 major and minor axis coordinates (X and Y, respectively). The PNe designated 
bv lMerrett et al.. I.200&) as "Stream" or "Stream?" are marked by red diamonds, and overlay the analogous stream in the simulation. 



test of the scenario in this paper, since it links the outbound 
NE shelf stars to the W shelf. 

The W shelf in our simulation is much less flattened 
than the NE shelf, and hence looks more like our idealized 
notion of a shell in Figure The majority of the stars lie 
within a cone 60° -90° wide, though many of the stars still lie 
on a curved concentration on one side of this cone. The in- 
bound and outbound stars are intermingled, and some of the 
inbound stars have been scattered far outside of the main 
conical region. The stream connecting the E and W stars 



passes very close ( < 1 kpc) to M31's center in this simula- 
tion, which helps account for the increased scatter in orbit 
directions. The main body of the shelf lies more distant than 
M31's center, and is viewed edge-on over much of its appar- 
ent outer boundary, which explains the relatively sharp and 
circular appearance in Figure |5| 

We note that there is yet another radial shell of debris 
in our model, consisting of particles even further forward 
than those in the W shelf. This shell is very faint, as it 
contains only about 6 x 10^ -^^0; which is only about a tenth 
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as much mass as the sheU making up the W shelf. This result 
is highly model-dependent, since this shell represents the 
extreme forward edge of the satellite debris. The shell may 
not exist at all, or may be much more massive than found 
in our simulation. The particles in this shell cover the entire 
E side of M31 fairly evenly, out to a radius of 18kpc. Deep 
kinematic surveys around the SE minor axis would seem to 
be the method most likely to detect this debris. 

4.3 Velocity pattern: expectations and relation to 
M31's potential 

Velocity measurements of the shell debris, as discussed by 
MK, are potentially the most powerful constraints on accept- 
able models of their makeup. This is partly because veloci- 
ties can be measured much more accurately than distances, 
and partly because the velocities are extremely coherent. 
Figure ITITk shows the location of the particles in the space 
of M31-centric radius versus M31-centric radial velocity. The 
particles occupy a cold stream, which spirals in the clock- 
wise direction as one goes further forward in the stream. 
The southern stream, NE shelf, and W shelf are easily seen 
as the first three groups of particles starting from the right. 
The fainter shelf just described on the eastern side is seen 
as the fourth group. The energy decreases along the stream, 
reducing the apocentric radii of successive radial loops. 

In qualitative accord with the predictions of MK dis- 
cussed above, we find that the velocities lie within a well- 
defined, roughly triangular region in the space of vios versus 
Rproj, as shown in Figure^] (see also the similar Figure^Ja). 
Over most of the region, the velocity distribution has two 
peaks at positive and negative velocity, which both approach 
II ~ (using velocity relative to M31) as Rproj approaches 
the shelf boundary. These two peaks are shown clearly in the 
histograms of the velocity distribution in Figure [TTI These 
two peaks generally have unequal strength, particularly at 
small radii in the W shelf where the negative-velocity peak 
is a fraction of its positive- velocity counterpart. The return- 
ing or positive-velocity "stream" in the NE shelf shows up 
quite clearly in this representation. As mentioned above, this 
component has so far eluded detection, probably due mostly 
to the areal coverage of the kinematic surveys. 

This velocity pattern will be easiest to recognize near 
the edges of the shelf, in the regions far away from M31's 
disk. We discussed evidence that portions of this pattern are 
visible in existing observations in § 3.3, including the appar- 
ent outbound "stream" in the NE shelf. The corresponding 
inbound stream is more easily seen in Figures ITni and lTTI than 
in Figure |3| both because the division at y = —0.3° splits 
the stream particles between the two panels of the latter 
figure, and because the line-of-sight velocity is a function of 
Rproj rather than X, and those two quantities are not well 
aligned in the main region occupied by the inbound stream. 

We fit the upper and lower envelopes in panels (b) and 
(c) of Figure lTUl with straight lines, using only the outer por- 
tions of each shelf. Using MK's approximation in Equation|H| 
we obtain estimated circular velocities of 291kms~^ in the 
NE shelf at 36 kpc from M31, and 245 kms"^ in the W shelf 
at 24 kpc. The actual circular velocities at these two radii in 
our potential model are 214 and 229 kpc, respectively. Thus, 
the MK prediction is not in precise quantitative accord with 
the simulation, but rather overestimates the enclosed mass 



by 15-85%. There are other disagreements with the MK 
description as well: the velocity slopes are not equal and 
opposite, and the tip velocity is nonzero. 

Using the simulation particles, we find that these dis- 
crepancies originate from several effects not considered by 
MK, including the following: the nonzero angular momen- 
tum in the debris; the significant gradient of energy along 
the stream; the continuing dynamical interaction with the 
progenitor, mainly in the NE shelf; perspective effects from 
the large size and small distance of the M31 system; and fi- 
nally, the finite range of radius needed to measure the slope. 
All of these effects change the velocities visibly in Figure [101 
and taken together change the measured slope significantly. 
They also break the above-mentioned symmetry between 
the caustics from particles on opposite sides of M31, though 
the particle distribution is so heavily weighted to one side in 
each case that this effect is barely visible. In both shells, the 
limited spatial coverage of the shell does not greatly modify 
the slope of the velocity envelope, although the lack of stars 
near the shell edge does deplete the region near the tip and 
make the slope measurement more challenging. 

Despite these complications, the general point of the 
MK paper is still valid: the velocity signature in shells can 
serve as a useful diagnostic of the gravitational potential, 
over a radial range 15-35 kpc in this case. Our gravitational 
potential is not necessarily the correct one. We expect the 
measurement biases to be similar for different potentials 
though, so to first order the ratio found here between pre- 
dicted and measured velocity gradients can be used to es- 
timate the true enclosed mass at the shelf radii. Measure- 
ments o f the H I rotation cu rve now probe out to ~ 35 kpc 
as well dCarignan et al ]|200d) . but there are significant con- 
cerns about these results due to the warping and irregular 
structure of the disk. Hence, the shelf velocities will help con- 
firm the H I measurements, and will also allow comparison 
of the potential in and out of the plane, thereby constraining 
the ellipticity of the dark matter distribution. 

4.4 Surface density profile as a diagnostic of the 
progenitor 

We now investigate what the surface density of the shelves 
signifies for the position and mass of the progenitor. First, 
we check the approximate relationship between the surface 
density and the mass fiow rate, using the simulated W shelf 
as a test case. Equation|S]is the appropriate form to use for 
this shelf since it is nearly tangent to the line of sight. In 
§ 3.2 we found the simulated surface density in the W shelf 
was Tjw ~ 7.2 X lO^MQkpc"^. For the other constants 
in Equation |S| we find rs ~ 25 kpc, Vc ~ 230kms~^, and 
dM/dt ~ 1.6AfQyr~^, using the rate of mass in the shell 
reaching turnaround as an estimate of the latter quantity. 
Converting the particles in this radial shell into a fully spher- 
ical shell by randomizing their orientations relative to M31, 
and computing the surface density in the W shelf region 
as before, we find a surface density of 1.5 x IO^A/q. This 
is about 5 times smaller than the actual simulated density, 
implying fe/fn = 5 in Equation |S1 which is the appropriate 
form to use for the W shelf since it is nearly tangent to the 
line of sight. The expected surface density of the W shelf 
is then ~ 6.3 x lO^AfQkpc"^, consistent with the actual 
surface density in the simulation within the limits of our 
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Figure 11. Histograms of the line-of-sight velocity of satellite debris particles at different projected radii. Left panel: Particles in the 
quadrant -X'j\,/3i < 0, Yj^i^i < containing the NE shelf. Curves show results for 0.5° < Rproj < 0.6° (black), 1.0° < Rproj < 1-1° 
(blue), 1.5° < Rproj < 1-6° (green), 2.0° < Rproj < 2.1° (red). Right panel: Same, for the quadrant Xm3i > 0, Ymzi > containing 
the W shelf. Curves show results for 0.5° < Rp 



< 0.6° (black), 1.0° 



(blue), 1.5° < Rproj < 1-6° (green). 



accuracy. Predicting the NE shelf surface density is poten- 
tially more complicated, since it may contain contributions 
from both the first and second pericentric passages of the 
progenitor satellite, so we will not make use of it here. 

By viewing successive snapshots of our A''-body sim- 
ulation, we can easily see that the surface density in the 
W shell increases with time, because the mass flow rate in 
the stream dM/dt is also increasing as the progenitor ap- 
proaches closer to this point in the orbit. This suggests that 
we can use the shelf surface brigh tness to c onstrain the pro- 
genitor mass and orbital phase. In lPaoer ij. we put limits on 
the mass and phase using properties of the southern stream 
alone. The basic idea was that given a progenitor phase, 
the progenitor mass sets both the overall normalization and 
the phase extent of the stream mass, and hence the fraction 
of the mass in the visible extent of the southern stream. 
This argument was limited in its power, though, because a 
larger progenitor phase could always be compensated for by 
a larger progenitor mass. With the detection of the W shelf, 
it seems that the progenitor phase can be limited from the 
other direction, effectively triangulating its position. 

We proceed similarly to §3.2 in iPaoer III We fit a se- 
quence of trajectories assuming different present-day pro- 
genitor phases Fp, and otherwise following the method of 
§3.1 above. For each fit, we calculate the time scali ng factor 
r at the W shelf and at Fields 4 and 8 oflMcConnac hie et alJ 



(|2003D within the southern stream, wh i ch spa n the same ra- 
dial range as that used bv llbata et al.l (120011) in calculating 
the stream's luminosity. We also calculate the progenitor's 
orbital energy Eq, its progenitor's radial period tro, its peri- 
centric radius Vp, the W shell radius Vs, and the circular 
velocity at rp and rs. 

For the distribution of orbital energy in the debris, 
x{E), we adopt a Gaussian form with the dispersion as 
given in Equation This form is a reasonably accurate ap- 



proximation to the simulations in iPaoer III and this paper, 
and an improvement on the top hat distribution used in 
Fardal et al. We estimate the energies at Fields 4 and 8 as 

1/3 1/3 

E4 = Eqt^ and Eg = Eqt^ , using the scalings above 
with k — —0.4. We take the total mass of the stream be- 
tween fields 4 and 8 as Mstream ~ 2 x 10* Mq, as discussed 
in § 13.21 We then estimate the mass of the progenitor satel- 
lite Ms from that of the stream by solving the following 
equation: 

Ms / x{E)dE^Ms tream • (9) 
J Eg 

This result is shown by the solid curve in Figure m\ 

We then estimate Ms independently from the surface 
density of the W shelf using Equation |S| We relate the mass 
flow rate at the edge of the shelf to the steadily increasing 
energy there using 

dM , ^ rT.^dE 
^-^^^(^^rfT' 

where the change in energy per unit time t at a flxed orbital 
phase is 

dE _ dE dr ^ E 
dt dr dt 3t 

Combining these results with Equation|n] we get an equation 
that can be solved for the progenitor satellite's mass AIs'. 

fn 



Ms^(rsn)x{E) 



2--^r shVc{r sh)^{r sh) 

JB 



(10) 



Here we take fe/ fn = 5, as found using the simulated de- 
bris. In Figure IT^ we provide four curves using shell surface 
densities of "Ew = 2, 4, 8, and 16 x lO^Af0kpc"^, which 
probably cover the allowed range (8 x 10^ is our most likely 
value, although this includes the smooth halo component of 
M31). 
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Figure 10. Panel a: M31-contric radial velocity versus M31- 
centric radius for the satellite debris particles in the simulation. 
The southern, trailing stream is made up of the long tail to large 
radius, while the debris further forward spirals clockwise in a 
series of radial loops. The NE shelf is made up of the particles in 
the largest complete loop, out to 40 kpc, while the W shelf is made 
up of the particles in the next loop. Panel b: Line-of-sight velocity 
versus projected radius, showing only particles in the quadrant 
^M3i < 0, Ym31 < containing the NE shelf. Particles genuinely 
in the radial shell corresponding to the NE shelf are marked with 
red dots. The black dots show particles belonging to other radial 
shells; these are mostly spillover from the southern stream, though 
some debris further forward than the W shelf is also present at 
small radii. Note that in this projection, the returning "stream" 
in the NE shelf (the concentration at the upper edge of the red 
particles) is more easily seen than in Figure |21 Panel c: Same 
as Panel b, but showing particles in the quadrant Xj(/3i > 0, 
Ym3i > containing the W shelf. Particles genuinely in the radial 
shell corresponding to the W shelf are marked with red dots. 
Other particles, mostly spillover from the NE shelf, are shown as 
black dots. 
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Figure 12. Estimates of the progenitor's mass at the time it dis- 
rupts to form the stream, as a function of the progenitor's current 
orbital phase. The orbital phase is defined in units of the progen- 
itor's radial orbital period, starting from the time of disruption. 
The derivations of the various estimates are discussed in the text. 
Rising solid curve: estimate from stellar mass in the stream, ac- 
counting for the mass outside the stream. Falling curves: estimates 
from the surface density in the W shelf, using surface densities of 
2, 4, 8, and 16 X 10^ Mq kpc~^ from bottom to top. Horizontal 
dashed line: estimate of stellar mass in the stream alone. Horizon- 
tal dot-dash ed line: estimate from the metallicity of stream stars 
jFont et al.k2006) . With a better-constrained W shelf surface den- 
sity, the intersection of the rising and falling curves can be used 
to constrain the progenitor mass and current orbital phase. 



The estimate of the progenitor mass and phase is then 
given by the point where the W shelf curve, for an observed 
T,w, crosses the solid curve from the stream mass. Both an- 
alytic estimates use a series of approximations and may be 
imprecise by perhaps a factor of two. Also, assuming that 
the stripped fraction or baryonic fraction in the progenitor 
is less than unity will raise both mass estimates together. 
The plot suggests a minimum progenitor mass of Ms > 1.5- 
5 X 10^ Mq , depending on the shelf surface density. These 
masses are consistent with our simulation and with the esti- 
mate of the progenitor mass from the metallicity { Fo nt et alJ 



1200(1) ■ shown as the upper horizontal line. On the other 
hand, these estimates are far above the integrated mass in 
the stream itself, shown as the lower horizontal line. 

The progenitor's orbital phase can already be limited to 
Fp > 1.1, according to the plot. The plot also shows that Fp 
can be determined within a small fraction of a radial period, 
once the surface density in the W shelf is known precisely. 

The constraints given here can be improved in the fu- 
ture by using surface density constraints all the way out to 
the tip of the southern stream, and by using an ensemble of 
A^-body simulations to avoid the need for analytic approxi- 
mations. Nevertheless, it should be clear from this discussion 
that the surface density distribution is an important clue to 
the orbital and physical parameters of the progenitor. 
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5 DISCUSSION 

We have performed a simple A^'-body simulation in the grav- 
itational potential of M31, and used this to show that an 
accreting satellite can explain the giant southern stream, 
the eastern shelf, and a similar shelf on the western side 
as a single continuous structure. Following the entire folded 
path traced by the stream and shelf particles, the struc- 
ture is > 230 kpc long, certainly among the largest coher- 
ent stellar structures known despite its modest total stel- 
lar mass of ~ 2 X 10^ Mq. The NE and W shelves in this 
model represent successive radial loops of the orbit broad- 
ened into fan shapes by passage through pericenter. The 
NE shelf lies closer to us than the disk of M31, and partially 
overlaps with it. This may have relevance to microlensing 
surve ys which aim to detect compact objects of stella r-type 
mass JCalchi Novati et aLllioollde Jong et alj|200'i . since 
the lensing event probability will be boosted relative to disk 
self-lensing by the gap in distance between the disk and NE 
shelf, and by the large relative velocities between the disk 
and shelf stars. 

The counter-rotating "stream" in the NE seen in both 
PNe and RGB stars in this view represents a continuation 
of the NE shelf onto the face of M31, while the apparently 
narrow velocity width of this stream is explained by a caus- 
tic in velocity space. The properties of the incoming satel- 
lite are reasonable for a galaxy of its initial stellar mass, 
as discussed in § 3.1. While we have not explored velocity 
distributions outside the shelf regions in detail, our impres- 
sion is that the model is satisfactory in terms of where it 
does not put debris, as well as where it does. O ur expe- 
rience with different initial conditions er IIT suggests 
that this would not necessarily be true for orbits that tried 
to associate the southern stream with other objects such as 
M32 or the Northern Spur. 

We experimented with low-resolution versions of our 
main run, testing different orbits and satellite properties. 
Our sense acquired during this exploration is that the vol- 
ume of parameter space that can match current observations 
is large enough not to need extreme fine-tuning. The time 
corresponding to "present day" within our run is uncertain 
by perhaps 100 Myr, and it is not clear whether the progen- 
itor need be disrupted or not. As surveys of M31's disk and 
halo regions continue, we expect that the requirements on 
the model will become much more stringent. 

We remind the reader that we have neglected several 
physical effects in computing this model that are potentially 
important. These include dynamical friction, the response 
of the modes in the disk, rotation of the accreting satellite, 
dark matter in and around its stellar component, and the 
hydrodynamic interaction of gas within the satellite with the 
disk and halo of M31. In addition, the observations do not 
yet cover some of the most interesting areas for this model, 
and our comparison to the existing observations is as yet 
rudimentary. We therefore regard the existing simulation as 
a useful illustration, rather than a definitive prediction of the 
debris structure. Work on the interaction of the satellite with 
live dynamical components of M31 is currently underway. 

Kinematic observations in the shelf regions are impor- 
tant for two reasons. First, they can help clarify the nature 
and extent of the debris related to the giant southern stream. 
The many testable predictions of our model include: 



• The "PNe stream" of iMerrett et alJ J200d) should be 
continuously connected to the NE shelf, since in our model 
it is merely a region where the kinematics in projection gives 
a roughly constant velocity, and not a separate structure. 

• In both shelves, the characteristic kinematic pattern of 
a shell should be evident in the Rproj-vios plane (Figure^J. 

• There should be stars with positive velocities returning 
to the center of M31 in the NE shelf, which connect the NE 
and W shelves. 

• There may be a yet weaker shelf on the E side, from 
debris further forward in the stream. 



Second, the velocity signature of the shells depend on the 
potential of M31. Many more spectroscopic fields will be 
required to allow a precise measurement of the potential. 
However, unlike the case of the coUimated southern stream, 
the line-of-sight distances need not be computed at all to 
probe the potential of M31, as the geometry itself selects 
the stars in the velocity caustic. This removes one of the 
main sources of uncertainty in constraining the mass model. 
Comparison of the potential derived from the shelves to that 
obtained from H I rotation curves can help test the deviation 
of M31's mass distribution from spherical symmetry. 

If confirmed as shells connected to the southern stream, 
the NE and W shelves would to our knowledge represent 
the first detection of a multiple-shell system in a late-type 
spiral galaxy. Shells are very common in elliptical galax- 
ies. 137 elliptica l gala xies wit h shells were catalog u ed by 
iMalin fc Carteil lll983l) . while iMerrifield fc KuiikenI lll99ai 
estimate that 10% of elliptical galaxies have shells. The rea- 
sons for this strong morphological bias have not been com- 
pletely explained at present. Spirals may be destroyed and 
ellipticals created in the mergers that produced the shells, 
although this would only be plausible for merging systems 
much more massive than in the case discussed here (see 
iHernauist fc SoergellQQ^ . for a discussion of shell formation 
in major mergers). Alternatively, shells may dissolve faster 
in the non-spherical potential of the disk due to differential 
precession of the orbits. The difficulty of distinguishing shells 
from spiral features in late- type galaxie s may also contribute 
to th e lower rate of observing shells iBarnes fc Hernauisd 
Il992h . If all shells around spiral galaxies were as faint as the 
ones in M31, they would have escaped detection, sin ce RGB 
star count maps like those of iFerguson et al.l ^2002^ can be 
made only at close range. M31's putative shell system would 
also be the first in which kinematic observations have been 
made, allowing us to compare to and eventually constrain 
mass models of the host galaxy. 

Theoretical studies of dark halos have established their 
mean mass profile to high accuracy in the absence of 
baryons. Their response to the baryonic component in terms 
of their mass profile, ellipticity, and alignment with the disk 
is more difficult to simulate, making observational checks 
important. However, these aspects are difficult to study in 
most galaxies. The abundance in M31 of both coherent trac- 
ers of the potential, such as the tidal debris discussed in this 
paper, and statistical tracers such as RGB stars, satellite 
galaxies, and globular clusters, suggests that this galaxy can 
be a Rosetta Stone for understanding the properties of dark 
halos and their interaction with the galaxies they host. 
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